#
# This is a Shiny web application. You can run the application by clicking
# the 'Run App' button above.
#
# Find out more about building applications with Shiny here:
#
#    http://shiny.rstudio.com/
#
#options(repos = BiocManager::repositories())
#options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
#BiocManager::install(c(
 # "limma","GEOquery","BiocGenerics","Biobase"
#))
library(shiny)
library(dplyr)
library(fs)
library(stringr)
suppressPackageStartupMessages(library(GEOquery)) 
#gset=AnnoProbe::geoChina('GSE11318')
#load("data/GSE11318_eSet.Rdata")
library(AnnoProbe)
library(limma)
library(DT)

# 定义本shiny的基本函数



#基因表达矩阵


get_express_matrix <- function(gset){
  probes_expr <- gset[[1]] %>%  exprs()
  probes_expr <- probes_expr+1  %>% log2()
}

# 临床表型数据

get_pho_data <- function(gset){
  pdata <- pData(gset)
}

# 获取基因注释信息及探针

get_gpl_probe <- function(probes_expr,gset){
  probe2gene <- gset@annotation %>% idmap()
  genes_expr <- filterEM(probes_expr,probe2gene)
}

# make group list 为差异基因分析做准备
# 在此之前，进行抽样，使得数据平衡
# 在group里面抽样使得数据平衡
sample_balanced_inside <- function(clinical_beataml_BM,types,ratio = 1){
  long_df <- clinical_beataml_BM %>% filter(group == types[1]) 
  short_df <- clinical_beataml_BM %>% filter(group == types[2]) 
  
  #sam <- sample(1:nrow(short_df),nrow(long_df))
  #short_df <- short_df[sam,]
  if (nrow(long_df) >= nrow(short_df)){
    sam <- sample(1:nrow(long_df),nrow(short_df)*ratio)
    long_df <- long_df[sam,]
  }else{
    sam <- sample(1:nrow(short_df),nrow(long_df)*ratio)
    short_df <- short_df[sam,]
  }
  df <- rbind(long_df,short_df)
}

get_group_samples <- function(gset,types){
  pdata <- pData(gset) %>% select(characteristics_ch1.6,geo_accession) %>% 
    filter(characteristics_ch1.6 != "Clinical info: Final microarray diagnosis: Unclassified DLBCL") %>% 
    mutate(group = case_when(
      characteristics_ch1.6 == "Clinical info: Final microarray diagnosis: PMBL" ~ "PMBL",
      characteristics_ch1.6 == "Clinical info: Final microarray diagnosis: GCB DLBCL" ~ "GCBDLBCL",
      characteristics_ch1.6 == "Clinical info: Final microarray diagnosis: ABC DLBCL" ~ "ABCDLBCL"
    ) ) %>% 
    filter(group %in% types) %>% 
    sample_balanced_inside(types,ratio = 1) %>% 
    rename(sample_ID = geo_accession) %>% 
    select(-characteristics_ch1.6) %>% 
    relocate(sample_ID) %>% 
    arrange(group)
}

get_group_design <- function(group_sample){
  group_list=factor(group_sample$group)
  table(group_list)
  design=model.matrix(~factor(group_list))
  return(design)
}

get_expression_sort <- function(gset,group_sample){
  probes_expr <- exprs(gset)
  probes_expr <- probes_expr+1  %>% log2()
  
  probe2gene <- gset@annotation %>% idmap()
  genes_expr <- filterEM(probes_expr,probe2gene) %>% select(all_of(group_sample$sample_ID))
  
  return(genes_expr)
  
}
# 进行deg分析

perform_DEG_analysis <- function(genes_expr,design){
  genes_expr %>% 
    lmFit(design) %>% 
    eBayes() %>% 
    topTable(coef=2,n=Inf)
}

# DEG workflow

DEG_workflow <- function(gset,types){
  
  group_sample <- gset %>% get_group_samples(types = types)
  design <- group_sample %>% get_group_design()
  DEG <- gset %>% 
    get_expression_sort(group_sample) %>% 
    perform_DEG_analysis(design)
  return(DEG)
}

## visualization

deg_volcano_graph <- function(DEG,id){
  #DEG <- DEG[[ids]]
  need_deg=data.frame(symbols=rownames(DEG), logFC=DEG$logFC, p=DEG$P.Value)
  deg_volcano(need_deg,id)
}

deg_heatmap_graph <- function(DEG,genes_expr,group_list,id){
  #DEG <- DEG[[ids]]
  #genes_expr <- genes_expr[[ids]]
  #group_list <- group_list[[ids]]
  deg_heatmap(DEG,genes_expr,group_list,id)
}

# make DEG tables

deg_tables <- function(need_deg,p_thred=0.05,logFC_thred=1){
  need_deg <- need_deg %>% 
    select(all_of(c("logFC","adj.P.Val"))) %>%
    rename(p = adj.P.Val) %>% 
    mutate(change = ifelse(p < p_thred & abs(logFC) > logFC_thred,
                                       ifelse(logFC > logFC_thred ,'UP','DOWN'),'NOT')) %>% 
    arrange(p) %>%
    na.omit()
}

# Define UI for application that draws a histogram
ui <- fluidPage(

    # Application title
    titlePanel("GSE11318 DEG Utils"),

    # Sidebar with a slider input for number of bins 
    sidebarLayout(
        sidebarPanel(
            sliderInput("seed",
                        "Set Seeds:",
                        min = 1,
                        max = 100,
                        value = 42),
            textInput("GSE_id",
                      "Tell me the GSE IDs:",
                      "GSE11318"),
            textInput("target_varibale",
                        "Select Your Group Varible: ",
                        "工事中，主不在乎"),
            textInput("magic_formular",
                      "Prompt Your Magic Formular",
                      "PMBL ~ ABCDLBCL"
                      ),
            textInput("gene_id",
                      "Write down Your dream of Genes",
                      "CCL5")
        ),

        # Show a plot of the generated distribution
        mainPanel(
          tabsetPanel(
            tabPanel(title = "DEGs表",
                     DTOutput("DEGs_table")
                     ),
            tabPanel(title = "火山图",
                     plotOutput("volcano"),
                     plotOutput("volcano2")
            ),
            tabPanel(title = "热图",
                     plotOutput("hetmaps")
            ),
            tabPanel(title = "兴趣基因检验", 
                     plotOutput("gene_diff"),
                     tableOutput("gene_table")
            )
          )
        )
    )
)

# Define server logic required to draw a histogram
server <- function(input, output) {
  
  gset <- reactive({
    withProgress(message = 'load datas ...',
                 detail = '别急，急也没有用...', value = 1, {
                   
    if (file_exists("data/GSE11318_eSet.Rdata")){
      load("data/GSE11318_eSet.Rdata")
    } else {
      gset = AnnoProbe::geoChina('GSE11318')
    }
    gset <- gset
    
  })
  })
  
  types <- reactive({
    set.seed(input$seed)
    words <- str_split(input$magic_formular,"~") %>% 
      unlist()
    for (i in 1:length(words)) {
      words[i] <- str_trim(words[i])
    }
    words <- words
  })
  
    DEG <- reactive({
      gset <- gset()
      set.seed(input$seed)
      DEG_workflow(gset[[1]],types = types())
    })
    
    group_sample <- reactive({
      gset <- gset()
      set.seed(input$seed)
      gset[[1]] %>% get_group_samples(types = types())
    })
    genes_expr <- reactive({
      gset <- gset()
      set.seed(input$seed)
      group_sample <- group_sample()
       gset[[1]] %>% 
        get_expression_sort(group_sample)
    })
    gene_list <- reactive({
      set.seed(input$seed)
      group_sample <- group_sample()
      gene_list <- factor(group_sample$group)
    })
    
    
    output$DEGs_table <- renderDT({
      deg_tables(DEG(),p_thred=0.05,logFC_thred=1)
    })
    output$volcano <- renderPlot({
      deg_volcano_graph(DEG(),2)
    })
    
    output$volcano2 <- renderPlot({
      deg_volcano_graph(DEG(),1)
    })
    
    output$hetmaps <- renderPlot({
      withProgress(message = '正在为您作画，可能有点久 ...',
                   detail = '别急，急也没有用...', value = 1, {
      #deg_heatmap_graph(DEG(),
       #                 genes_expr(),
        #                gene_list(),
         #               10)
      library(pheatmap)
      deg <- DEG()
      genes_expr <- genes_expr()
      group_list <- gene_list()
      x = deg[,1]
      topn <- 20 
      names(x)=rownames(deg)
      cg=c(names(head(sort(x),topn)),names(tail(sort(x),topn)))
      n=t(scale(t(genes_expr[cg,])))
      n[n>2]=2
      n[n< -2]= -2
      n[1:4,1:4]
      ac=data.frame(group_list=group_list)
      rownames(ac)=colnames(n)
      p <- pheatmap(n,show_colnames =F,show_rownames = T,annotation_col=ac)
      library(ggplotify)
      as.ggplot(p) + theme(plot.margin = margin(l = 5, r = 10))
                   })
    })
    
    output$gene_diff <- renderPlot({
      check_diff_genes(input$gene_id,
                       genes_expr(),
                       gene_list())
    })
    
    output$gene_table <- renderTable({
      DEGtable <- DEG()
      DEGtable[input$gene_id,]
    })
}

# Run the application 
shinyApp(ui = ui, server = server)
